
plot_episodes <- function(years = c(1900, 2021),
                          country = c("Sweden"),
                          data = episode_data_wo_CI) {
  
  eps <- episode_data_wo_CI
  
  eps <- eps %>%
    dplyr::group_by(country_id) %>%
    mutate(increase_episode = ifelse(lead(increase_episode == 1), 1, increase_episode), 
           decline_episode = ifelse(lead(decline_episode == 1), 1, decline_episode))
  
  stopifnot(is.numeric(years), length(years) == 2, years[2] > years[1])
  
  stopifnot(is.character(country))
  
  if(length(country) > 1)
    stop("Error: More than one country selected")
  
  if(length(country) == 0)
    stop("Error: No country selected")
  
  if(!country %in% data$country_name)
    stop("Error: Country not found")
  
  if(max(years) < min(eps %>% filter(country_name==country) %>% pull(year)) | max(years)>max(eps %>% filter(country_name==country) %>% pull(year)))
    stop("Error: Data not available for time range")
  
  year <- country_name <- increase_episode <- decline_episode <- overlap_eps <- v2xca_academ <-
    ep_type <- episode <- vdem <- increase_episode_start_year <-increase_episode_end_year <-
    decline_episode_start_year <- decline_episode_end_year <- decline_episode_id <- increase_episode_id <- countries <- NULL
  
  eps_year <- eps %>%
    dplyr::filter(country_name == country, dplyr::between(year, min(years), max(years))) %>%
    dplyr::filter(increase_episode == 1 | decline_episode == 1) 
  
  if(nrow(eps_year)>1){
    eps_year <- eps_year %>% 
      dplyr::mutate(overlap_eps = ifelse(!is.na(decline_episode_id) & !is.na(increase_episode_id), "overlaps", NA)) %>% 
      tidyr::pivot_longer(cols = c(decline_episode_id, increase_episode_id, overlap_eps), names_to = "ep_type", values_to = "episode") %>%
      dplyr::select(country_name, year, v2xca_academ, ep_type, episode,
                    decline_episode_start_year, decline_episode_end_year, 
                    increase_episode_start_year, increase_episode_end_year) %>%
      dplyr::filter((ep_type == "increase_episode_id") |
                      (ep_type == "decline_episode_id") |
                      ep_type == "overlaps") %>%
      drop_na(episode) %>%
      group_by(year) %>%
      mutate(overlap_eps = n(),
             episode_id = ifelse(ep_type == "decline_episode_id", paste0("Decline: ", decline_episode_start_year, "-", decline_episode_end_year), episode),
             episode_id = ifelse(ep_type == "increase_episode_id", paste0("Growth: ", increase_episode_start_year, "-", increase_episode_end_year), episode_id)) %>%
      ungroup()
    
    v2xca_academ <- eps %>%
      filter(country_name == country, between(year, min(years), max(years))) %>%
      ungroup() %>%
      select(year, v2xca_academ)
    
    if(max(eps_year$overlap_eps) > 1) {
      print("Warning: Some episodes overlap!")
    }
    
    p <-   ggplot2::ggplot() +
      geom_line(data = eps_year, aes(group = episode_id, color = episode_id, linetype = ep_type,x = year, y = v2xca_academ)) +
      geom_line(data = v2xca_academ, aes(x = year, y = v2xca_academ), alpha = 0.35) +
      scale_colour_grey(breaks = levels(factor(eps_year$episode_id[eps_year$episode_id!="overlaps"])),
                        name = "Episode", start = 0.01, end = 0.01) +
      scale_linetype_manual(name = "Episode type", breaks = c("decline_episode_id", "increase_episode_id", "overlaps"),
                            labels = c("Decline", "Growth", "Overlap"),
                            values = c("dashed", "dotted", "solid")) +
      scale_x_continuous(breaks = seq(round(min(years) / 10) * 10, round(max(years) / 10) * 10, 10)) +
      xlab("Year") +  ylab("Academic Freedom Index") + ylim(0,1) +
      theme_bw() +
      guides(color = guide_legend(override.aes = list(size = 0))) +
      ggtitle(sprintf("%s", country))
    
    if (isTRUE(length(which(eps_year$ep_type == "increase_episode_id")) > 0)){
      
      if (any(eps_year$year%in%c(eps_year$increase_episode_start_year))) {
        p <- p +  geom_point(data = eps_year, aes(x = year, y = ifelse(year == increase_episode_start_year, v2xca_academ, NA)), shape = 2, alpha = 0.75) 
        
      } else {
        p
      }
      
      if (any(eps_year$year%in%c(eps_year$increase_episode_end_year))) {
        p <- p +geom_point(data = eps_year, aes(x = year, y = ifelse(year == increase_episode_end_year, v2xca_academ, NA)), shape = 17, alpha = 0.75)
      } else {
        p
      }
    }
    
    if (isTRUE(length(which(eps_year$ep_type == "decline_episode_id")) > 0)) {
      
      if (any(eps_year$year%in%c(eps_year$decline_episode_start_year))){
        p <- p +  geom_point(data = eps_year, aes(x = year, y = ifelse(year == decline_episode_start_year, v2xca_academ, NA)), shape = 1, alpha = 0.75) 
      } else {
        p
      }
      if (any(eps_year$year%in%c(eps_year$decline_episode_end_year))){
        p<- p+ geom_point(data = eps_year, aes(x = year, y = ifelse(year == decline_episode_end_year, v2xca_academ, NA)), shape = 16, alpha = 0.75)
      } else {
        p
      }
    }
    p
    
    
  } else {
    print("No episodes during selected period.")
    
    polyarchy <- eps %>%
      filter(country_name == country, between(year, min(years), max(years))) %>%
      ungroup() %>%
      select(year, v2xca_academ)
    
    p <-ggplot2::ggplot() +
      geom_line(data = v2xca_academ, aes(x = as.numeric(year), y = v2xca_academ), alpha = 0.35) +
      scale_x_continuous(breaks = seq(round(min(years) / 10) * 10, round(max(years) / 10) * 10, 10)) +
      xlab("Year") +  ylab("Academic Freedom Index") + ylim(0,1) +
      theme_bw() +
      ggtitle(sprintf("%s", country))
    
    p
    
  }
}